Feature/higher order mesh particle interface #265
Conversation
…to reorder how it sees fit
…mentation - needs testing
…he original reference position
will-saunders-ukaea
added
Awaiting Review
|
@oparry-ukaea to review bary and geometry transport please all confirm that the tests pass. |
|
All tests pass using gcc 13.3 and adaptivecpp 24.02.0 With oneapi 2025.0.0 Are they passing for you @oparry-ukaea @seimtpow ? |
Are you getting 'KILLED BY SIGNAL: 4 (Illegal instruction)'? |
Thread 9 "unitTests" received signal SIGILL, Illegal instruction. (from gdb) |
|
Sounds like the opencl bug Will found before Christmas (see the standups slack) |
Ah I missed that - second one works for me. Thanks! |
| @@ -0,0 +1,4082 @@ | |||
| <?xml version="1.0" encoding="utf-8" ?> | |||
There was a problem hiding this comment.
Maybe just commit the compressed version of the mesh xml?
| * @param[in] coord The evauation point in the dimension of interest. | ||
| * @param[in] z_values A length num_phys array containing the quadrature | ||
| * points. | ||
| * @param[in] z_values A length num_phys array containing the quadrature |
There was a problem hiding this comment.
| * @param[in] z_values A length num_phys array containing the quadrature | |
| * @param[in] bw_values A length num_phys array containing the quadrature |
| * @param[in] coord The evauation points, size N, in the dimension of interest. | ||
| * @param[in] z_values A length num_phys array containing the quadrature | ||
| * points. | ||
| * @param[in] z_values A length num_phys array containing the quadrature |
There was a problem hiding this comment.
| * @param[in] z_values A length num_phys array containing the quadrature | |
| * @param[in] bw_values A length num_phys array containing the quadrature |
| static_assert(!Newton::local_memory_required< | ||
| Newton::MappingQuadLinear2DEmbed3D>::required, | ||
| "Did not expect local memory to be required for this Newton " | ||
| "implemenation"); |
There was a problem hiding this comment.
| "implemenation"); | |
| "implementation"); |
| @@ -79,22 +104,13 @@ class PackedGeom2D { | |||
|
|
|||
| // curve of the edge | |||
| auto curve = seg_geom->GetCurve(); | |||
| ASSERTL0(curve == nullptr, "Not implemented for curved edges"); | |||
| // A curve with n_points = -1 will be a taken as non-existant. | |||
There was a problem hiding this comment.
| // A curve with n_points = -1 will be a taken as non-existant. | |
| // A curve with n_points = -1 will be a taken as non-existent. |
| auto xmap = geom->GetXmap(); | ||
| const bool linear = (xmap->GetBasisNumModes(0) == 2) && | ||
| (xmap->GetBasisNumModes(1) == 2) && | ||
| (xmap->GetBasisNumModes(2) == 2); | ||
| if (linear) { |
There was a problem hiding this comment.
| auto xmap = geom->GetXmap(); | |
| const bool linear = (xmap->GetBasisNumModes(0) == 2) && | |
| (xmap->GetBasisNumModes(1) == 2) && | |
| (xmap->GetBasisNumModes(2) == 2); | |
| if (linear) { | |
| if (geometry_is_linear(geom) { |
Use geometry_is_linear from utility_geometry.hpp?
| @@ -341,7 +415,8 @@ TEST(ParticleGeometryInterface, CoordinateMapping3D) { | |||
| xi1[dimx] = dh_xi.h_buffer.ptr[dimx]; | |||
| eta0[dimx] = dh_eta.h_buffer.ptr[dimx]; | |||
| } | |||
| geom->GetXmap()->LocCollapsedToLocCoord(eta0, xi0); | |||
| // geom->GetXmap()->LocCollapsedToLocCoord(eta0, xi0); | |||
There was a problem hiding this comment.
| // geom->GetXmap()->LocCollapsedToLocCoord(eta0, xi0); |
| @@ -367,7 +442,9 @@ TEST(ParticleGeometryInterface, CoordinateMapping3D) { | |||
| eta1[dimx] = dh_eta.h_buffer.ptr[dimx]; | |||
| xi0[dimx] = dh_xi.h_buffer.ptr[dimx]; | |||
| } | |||
| geom->GetXmap()->LocCoordToLocCollapsed(xi0, eta0); | |||
| // geom->GetXmap()->LocCoordToLocCollapsed(xi0, eta0); | |||
There was a problem hiding this comment.
| // geom->GetXmap()->LocCoordToLocCollapsed(xi0, eta0); |
| lambda_test_wrapper(geom, geom_test); | ||
| } else { | ||
| // Unknown shape type | ||
| ASSERT_TRUE(false); |
There was a problem hiding this comment.
| ASSERT_TRUE(false); | |
| FAIL() << "Unknown shape type."; |
?
| GeometryInterface::Triangle geom_test{}; | ||
| lambda_test_wrapper(geom, geom_test); | ||
| } else { | ||
| ASSERT_TRUE(false); |
There was a problem hiding this comment.
| ASSERT_TRUE(false); | |
| FAIL() << "Expected only triangles or quads but found a different shape type."; |
?
|
All tests passing for me with GCC 11.3.0 / AdaptiveCpp 24.06.0 With OneAPI 2024.1.0/ DPC++ 2022.1.0 I get intermittent failures of It's always very marginal, e.g. Could/should the tolerance be relaxed a bit? |
| x0, x1, x2, num_phys0, num_phys1, num_phys2, physvalsv.data(), | ||
| div_space.data(), z0v.data(), z1v.data(), z2v.data(), bw0v.data(), | ||
| bw1v.data(), bw2v.data(), stride); | ||
| EXPECT_NEAR(to_test, to_test_stride, 1.0e-14); |
There was a problem hiding this comment.
This tolerance is a little too low for the some of the tests to pass for me with oneapi (1.0e-13 would be fine). Could it be set using the tol argument rather than hard-coded here?
| get_all_elements_3d(particle_mesh_interface->graph, geoms); | ||
| for (auto gx : geoms) { | ||
| std::pair<int, std::shared_ptr<Nektar::SpatialDomains::Geometry3D>> | ||
| tmp = {gx.first, gx.second}; |
There was a problem hiding this comment.
Why is tmp created instead of just using gx?
| this->num_bytes_local_memory = | ||
| std::max(this->num_bytes_local_memory, s); | ||
| }; | ||
| for (auto &geom : geoms_local) { |
There was a problem hiding this comment.
A structured binding might be good here
| } | ||
|
|
||
| inline void free_data_v(void *data_host) { | ||
| Generic3D::DataHost *h_data = static_cast<Generic3D::DataHost *>(data_host); |
There was a problem hiding this comment.
Does the DataHost that was created earlier at *data_host with the placement new get destroyed properly? This function frees the data in the buffer but I'm wondering if the memory occupied by the DataHost itself gets freed.
| std::unique_ptr<BufferDeviceHost<char>> dh_data; | ||
| std::unique_ptr<BufferDeviceHost<std::byte>> dh_data; | ||
| /// The data required to perform newton iterations for each geom on the host. | ||
| std::unique_ptr<BufferHost<std::byte>> h_data; |
There was a problem hiding this comment.
Maybe missing something, but I don't like the buffer of bytes/void pointer - doesn't seem necessary as everything is known at compile time.
For example each NEWTON_TYPE class could add e.g.
using data_type = DataDevice;
and then we could get that via typename NEWTON_TYPE::data_type and just make dh_data a pointer to the correct type - it could work something like this
template <class T> struct A {
int do_something(void *data) {
auto &under = static_cast<T &>(*this);
return under.do_something_v(static_cast<typename T::data_type *>(data));
}
void write(void *data) {
auto &under = static_cast<T &>(*this);
under.write_v(static_cast<typename T::data_type *>(data));
}
};
struct B : A<B> {
using data_type = double;
int do_something_v(data_type *data) { return (int)*data; }
void write_v(data_type *data) { *data = 1.0; }
};
struct C : A<C> {
using data_type = int;
int do_something_v(data_type *data) { return *data + 1; }
void write_v(data_type *data) { *data = 1; }
};
template <typename K> struct Xthing {
using KType = typename K::data_type;
KType *ptr;
Xthing() { ptr = new KType; }
~Xthing() { delete ptr; }
};
Also not sure the CRTP is needed AFAICS you should be able to just use static functions if there is no state in the derived classes.
Description
Type of change
Testing